# Code to generate list of barangays included in the matched subset

library(Matching)
library(AER)
library(ebal)
library(xtable)

## Load data
#load("bgy_geovars_merge_updated.RData")
#bgy <- read.csv("bgy_merged.csv", stringsAsFactors = F)


## Matching covariates
covars <- c("GEOCODE2", "cadt_bin_update2010",
            "logpop_pretrend",
            "ip_2000", "ethfrac_2000", "area_sqkm", "elev_mean", "elev_sd",
            "slope_mean", 
            "soil_index", "log_ncip_dist", "log_coast_dist", "log_road_dist", "religion_catholic_2000",
            "strong_mat_pretrend",
            "own_lot_pretrend",
            "elem_grad_pretrend",
            "leg_pretrend",
            "have_street_pattern_2000", "have_highway_2000", "have_church_2000",
            "have_market_2000", "have_elem_2000", "have_bhs_2000", "have_water_sys_2000", 
            "birth_reg_2000")

## Subset to eligible universe
bgy_universe <- bgy[which(bgy$universe_ncip_updated == 1 & bgy$reg_name != "CAR"), ]
bgy_match_data <- bgy_universe[, covars]

bgy_match_data_nona <- na.omit(bgy_match_data)

exclude <- c("GEOCODE2", "cadt_bin_update2010")
include <- covars[-which(covars %in% exclude)]
X <- bgy_match_data_nona[, include]

my_formula <- formula(paste("cadt_bin_update2010 ~", paste0(include, collapse = "+"), sep = " "))

set.seed(02139)
# Propensity score matching - logit
prop.mod.probit <- glm(my_formula,
                      data = bgy_match_data_nona, family = binomial(link = "probit"))
match_out_prop_probit <- Match(Tr = bgy_match_data_nona$cadt_bin_update2010, 
                              X = prop.mod.probit$fitted.values,
                              replace = F,
                              M = 1)
mb_probit <- MatchBalance(my_formula,
                         data = bgy_match_data_nona, match.out = match_out_prop_probit, print.level = 0)


matched_samp_prop <- c(match_out_prop_probit$index.treated, match_out_prop_probit$index.control)
bgy_match_data_prop <- bgy_match_data_nona[matched_samp_prop, ]
matched_samp_bgys <- bgy_match_data_prop$GEOCODE2

# Export
save(matched_samp_bgys, file = "matched_samp_bgys.RData")


